cg_regions <- read_rds(here::here("data/derived_data/cg_regions.rds"))
animal_regions <- read_rds(here::here("data/derived_data/animal_regions.rds"))
#prism_zip <- read_rds(here::here("data/derived_data/prism_zip.rds"))
cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
  filter(trait == "ww" & var == "cg_sol") %>%
  filter(n_animals > 4) %>%
  group_by(year) %>%
  mutate(year_sd = sd(value),
         year_mean = mean(value),
         keep = case_when(
           value > year_mean + year_sd ~ "YES",
           value < year_mean - year_sd ~ "YES",
           TRUE ~ "NO"
         )) %>% 
  ungroup() %>% 
  filter(keep == "YES") %>% 
#  sample_frac(0.2) %>%
      ggplot(aes(
        x = lng,
        y = lat,
        color = "goldenrod",
        size = n_animals
      )) +
      geom_point(
        alpha = 0.3,
        #size = 0.5
        ) +
      scale_size(
        range = c(1,4)
      ) +
      usa +
      scale_x_continuous(expand = c(0, 0)) +
      scale_y_continuous(expand = c(0, 0)) +
      coord_map("albers", lat0 = 39, lat1 = 45) +
      cowplot::theme_map() +
      labs(
        x = NULL,
        y = NULL,
        color = NULL,
        title = str_wrap( "Weaning weight CG solutions, low to high", width = 55)
      ) +
      #Set the "anchoring point" of the legend (bottom-left is 0,0; top-right is 1,1)
      #Put bottom-left corner of legend box in bottom-left corner of graph
      theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        # Padding around the legend
        # top, right, bottom, left
        # legend.box.margin = margin(b = 0, r = 0.2, unit = "cm"),
        # Padding around the plot
        plot.margin = margin(
          t = 0.7,
          r = 0,
          b = 0.7,
          l = 1,
          unit = "cm"
        ),
        plot.title = element_text(
          # size = 56,
          # family =  "lato",
          size = 22,
          vjust = 6,
          face = "italic"
        )
      ) +
      # https://stackoverflow.com/questions/32656553/plot-legend-below-the-graphs-and-legend-title-above-the-legend-in-ggplot2
      guides(
        color = FALSE,
        size = FALSE
      ) +
  gganimate::transition_reveal(value)
## 
Frame 1 (1%)
Frame 2 (2%)
Frame 3 (3%)
Frame 4 (4%)
Frame 5 (5%)
Frame 6 (6%)
Frame 7 (7%)
Frame 8 (8%)
Frame 9 (9%)
Frame 10 (10%)
Frame 11 (11%)
Frame 12 (12%)
Frame 13 (13%)
Frame 14 (14%)
Frame 15 (15%)
Frame 16 (16%)
Frame 17 (17%)
Frame 18 (18%)
Frame 19 (19%)
Frame 20 (20%)
Frame 21 (21%)
Frame 22 (22%)
Frame 23 (23%)
Frame 24 (24%)
Frame 25 (25%)
Frame 26 (26%)
Frame 27 (27%)
Frame 28 (28%)
Frame 29 (29%)
Frame 30 (30%)
Frame 31 (31%)
Frame 32 (32%)
Frame 33 (33%)
Frame 34 (34%)
Frame 35 (35%)
Frame 36 (36%)
Frame 37 (37%)
Frame 38 (38%)
Frame 39 (39%)
Frame 40 (40%)
Frame 41 (41%)
Frame 42 (42%)
Frame 43 (43%)
Frame 44 (44%)
Frame 45 (45%)
Frame 46 (46%)
Frame 47 (47%)
Frame 48 (48%)
Frame 49 (49%)
Frame 50 (50%)
Frame 51 (51%)
Frame 52 (52%)
Frame 53 (53%)
Frame 54 (54%)
Frame 55 (55%)
Frame 56 (56%)
Frame 57 (57%)
Frame 58 (58%)
Frame 59 (59%)
Frame 60 (60%)
Frame 61 (61%)
Frame 62 (62%)
Frame 63 (63%)
Frame 64 (64%)
Frame 65 (65%)
Frame 66 (66%)
Frame 67 (67%)
Frame 68 (68%)
Frame 69 (69%)
Frame 70 (70%)
Frame 71 (71%)
Frame 72 (72%)
Frame 73 (73%)
Frame 74 (74%)
Frame 75 (75%)
Frame 76 (76%)
Frame 77 (77%)
Frame 78 (78%)
Frame 79 (79%)
Frame 80 (80%)
Frame 81 (81%)
Frame 82 (82%)
Frame 83 (83%)
Frame 84 (84%)
Frame 85 (85%)
Frame 86 (86%)
Frame 87 (87%)
Frame 88 (88%)
Frame 89 (89%)
Frame 90 (90%)
Frame 91 (91%)
Frame 92 (92%)
Frame 93 (93%)
Frame 94 (94%)
Frame 95 (95%)
Frame 96 (96%)
Frame 97 (97%)
Frame 98 (98%)
Frame 99 (99%)
Frame 100 (100%)
## Finalizing encoding... done!

# anim_save(filename = here::here("figures/ww_solutions_appear.gif"), nframes = 48, fps = 2)

Summary: year-to-year CG solutions by region

  • “Your EPDs are not affected by weather/climate/environment”

Weaning weight

Contemporary group solutions

cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_line(
    effect_var = "cg_sol",
    trait_var = "ww",
    y_lab = "Mean CG solution",
    plot_title = "Weaning weight contemporary group solutions reflect year-to-year environmental trends"
    ) +
  lims(y = c(440, 580))

#ggsave(here::here("figures/ww_cg.solutions_line_year.png"), width = 10, height = 5.4)

Breeding value solutions

animal_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_line(
    effect_var = "bv_sol",
    trait_var = "ww",
    y_lab = "Mean breeding value",
    plot_title = "Weaning weight genetic trend: 1973-2019"
    )

#ggsave(here::here("figures/ww_bv.solutions_line_year.png"), width = 10, height = 5.4)

Maternal effect solutions

animal_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_line(
    effect_var = "mat_sol",
    trait_var = "ww",
    y_lab = "Mean maternal effect solution",
    plot_title = "Milk genetic trend: 1973-2019"
    )

#ggsave(here::here("figures/ww_mat.solutions_line_year.png"), width = 10, height = 5.4)

Post-weaning gain

(Removed data from 1972: weird outlier, only one region with data)

Contemporary group solutions

cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_line(
    effect_var = "cg_sol",
    trait_var = "pwg",
    y_lab = "Mean CG solution",
    plot_title = "Post-weaning gain contemporary group solutions reflect year-to-year environmental trends"
    ) 

#ggsave(here::here("figures/pwg_cg.solutions_line_year.png"), width = 10, height = 5.4)

Breeding value solutions

animal_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_line(
    effect_var = "bv_sol",
    trait_var = "pwg",
    y_lab = "Mean breeing value",
    plot_title = "Post-weaning gain genetic trend: 1973-2019"
    ) 

#ggsave(here::here("figures/pwg_bv.solutions_line_year.png"), width = 10, height = 5.4)
  • What’s up with 1972 region 5??
cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
      filter(trait == "pwg" & var == "cg_sol") %>%
      mutate(yr = lubridate::year(weigh_date),
             region = as.character(region)) %>%
      group_by(region, yr) %>%
  summarise(mean = mean(value)) %>% 
      bind_rows(
        cg_regions %>%
          filter(trait == "pwg" & var == "cg_sol") %>%
          mutate(yr = lubridate::year(weigh_date)) %>%
          group_by(yr) %>%
          summarise(mean = mean(value)) %>%
          mutate(region = "All regions")
      ) %>% 
  arrange(yr)

Year-to-year variance

  • Do stressful years result in a higher variance in CG solutions? I.e., more likely some producers will change management to compensate
    • Doesn’t appear so

Weaning weight

cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_var_line(
    effect_var = "cg_sol", 
    trait_var = "ww",
    plot_title = "The variance of weaning weight CG solutions changes little across regions or time"
  )

#ggsave(here::here("figures/ww_cg.solutions_variance.png"), width = 10, height = 12)

Post-weaning gain

cg_regions %>%
  filter(!region %in% c(4, 6)) %>%
  solutions_var_line(
    effect_var = "cg_sol", 
    trait_var = "pwg",
    plot_title = "The variance of post-weaning gain CG solutions changes little across regions or time"
  )

#ggsave(here::here("figures/pwg_cg.solutions_variance.png"), width = 10, height = 12)

Contemporary group solutions “over & under” the mean by year

  • Does need to be on a national/year basis? Or on a “overall mean for the region” basis (I.e., "is this year different than what’s normal for my region)?
  • Also maybe: plot only 1 +/- standard deviation? Take out points around the mean?

Weaning weight

Against regional mean

p1 <-
  cg_regions %>% 
  filter(!region %in% c(4, 6)) %>%
  filter(trait == "ww" & var == "cg_sol") %>%
  filter(n_animals > 4) %>% 
  mutate(yr = lubridate::year(weigh_date),
         yr = as.integer(yr)) %>%
  group_by(region) %>% 
  mutate(comp = mean(value)) %>% 
  ungroup() %>% 
  group_by(zip, yr) %>% 
  mutate(sz = sum(n_animals),
         zip_sol = mean(value)) %>% 
  ungroup() %>% 
  mutate(
         overunder = 
           case_when(
             zip_sol > comp ~ "Above mean",
             zip_sol < comp ~ "Below mean"
           ),
         howmuch = zip_sol - comp,
         howmuch = abs(howmuch)
           ) %>% 
  select(trait, var, yr, region, comp, zip, sz, overunder, howmuch, lng, lat) %>% 
  distinct() %>% 
  filter(howmuch < 300) %>% 
  solutions_map(
    effect_var = "cg_sol",
    trait_var = "ww",
    plot_title = "test",
    color_var = howmuch,
    size_var = sz,
    size_range = c(1, 4)
      ) +
  facet_wrap(~ overunder) +
  theme(legend.position = "bottom",
        strip.text.x = element_text(size = 16)) +
  gganimate::transition_time(yr) +
  labs(title = "Weaning weigh CG solutions: {frame_time}")
anim_save(filename = here::here("figures/ww_overunder_byregion.gif"), animation = plot, nframes = 48, fps = 1)

Against national/yearly mean

p2 <- 
  cg_regions %>% 
  filter(!region %in% c(4, 6)) %>%
  filter(trait == "ww" & var == "cg_sol") %>%
  mutate(yr = lubridate::year(weigh_date),
         yr = as.integer(yr)) %>%
  group_by(yr) %>% 
  mutate(comp = mean(solution)) %>% 
  ungroup() %>% 
  group_by(zip, yr) %>% 
  mutate(sz = sum(n_animals),
         zip_sol = mean(solution)) %>% 
  ungroup() %>% 
  mutate(
         overunder = 
           case_when(
             zip_sol > comp ~ "Above mean",
             zip_sol < comp ~ "Below mean"
           ),
         howmuch = zip_sol - comp,
         howmuch = abs(howmuch)
           ) %>% 
  select(trait, var, yr, region, comp, zip, sz, overunder, howmuch, lng, lat) %>% 
  distinct() %>% 
  filter(howmuch < 300) %>% 
  solutions_map(
    effect_var = "cg_sol",
    trait_var = "ww",
    plot_title = "test",
    color_var = howmuch,
    size_var = sz,
    size_range = c(1, 4)
    ) +
  facet_wrap(~ overunder) +
  theme(legend.position = "bottom",
        strip.text.x = element_text(size = 16)) +
  gganimate::transition_time(yr) +
  labs(title = "Weaning weigh CG solutions: {frame_time}")
#anim_save(filename = here::here("figures/ww_overunder_byyear.gif"), animation = plot_byyear, nframes = 48, fps = 1)

Weather

  • I think yearly temp and rainfall might not actually be as straightforward as I thought originally
    • I.e., a higher yearly temp could also mean a warm winter –> higher forage availability
      • See also: not taking into account differences in calving season across trends
    • Maybe eventually pull USDA drought data?
cg_regions %>%
  filter(var == "cg_sol" & trait == "ww") %>%
  mutate(year = lubridate::year(weigh_date)) %>%
  left_join(prism_zip %>%
              select(-lat,-lng) %>% 
              mutate(year = year - 1)) %>%
  group_by(year) %>%
  mutate(
    yr_sol = mean(value),
    yr_temp = mean(tmean),
    yr_rain = mean(ppt),
    sd_sol = sd(value),
    sd_temp = sd(tmean),
    sd_rain = sd(ppt),   
  ) %>%
  ungroup() %>%
  filter(region %in% c(1, 2, 5)) %>% 
  group_by(region, year) %>%
  summarise(
    dev_sol = (mean(value) - unique(yr_sol))/unique(sd_sol),
    dev_temp = (mean(tmean) - unique(yr_temp))/unique(sd_temp),
    dev_rain = (mean(ppt) - unique(yr_rain))/unique(sd_rain)
  ) %>%
  ungroup() %>%
  reshape2::melt(id = c("region", "year")) %>% 
  ggplot(
    aes(
      x = year,
      y = value,
      linetype = variable,
      color = forcats::as_factor(region)
    )
  ) +
  geom_line(size = 1.25) +
  scale_linetype_manual(
    values = c(
      "dev_sol" = "solid",
      "dev_temp" = "dotted",
      "dev_rain" = "twodash"
      ),
    labels = c(
      "dev_sol" = "CG solution",
      "dev_temp" = "Temperature",
      "dev_rain" = "Rain"
      )
    ) +
  scale_color_manual(
    values = c(
      "1" = "tomato2",
      "2" = "darkslategray4",
      "3" = "springgreen3",
      "4" = "brown",
      "5" = "goldenrod1",
      "6" = "gray50",
      "7" = "deeppink3",
      "8" = "gray17",
      "9" = "slateblue2"
    ),
    labels = c(
      "1" = "1: Desert",
      "2" = "2: Southeast",
      "3" = "3: High Plains",
      "4" = "4: Rainforest",
      "5" = "5: Arid Prairie",
      "6" = "6: Cold Desert",
      "7" = "7: Forested Mountains",
      "8" = "8: Fescue Belt",
      "9" = "9: Upper Midwest & Northeast"
    )
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(
      size = 22,
      face = "italic",
      margin = margin(
        t = 0,
        r = 0,
        b = 13,
        l = 0
      )
    ),
    axis.title = element_text(
      size = 16),
    axis.title.y = element_text(margin = margin(
      t = 0,
      r = 13,
      b = 0,
      l = 0
    )),
    axis.title.x = element_text(margin = margin(
      t = 13,
      r = 0,
      b = 0,
      l = 0
    )),
    axis.text = element_text(
      size = 14),
    legend.text = element_text(
      size = 14)
  ) +
  labs(x = NULL,
       y = "Standardized value",
       color = NULL,
       linetype = NULL,
       title = str_wrap("Environmental variables are not directly related to weaning weight CG solutions", width = 55)
       ) +
  facet_wrap(~ region, 
             nrow = 3) 
cg_regions %>%
  filter(var == "cg_sol" & trait == "ww") %>%
  mutate(year = lubridate::year(weigh_date),
         ) %>%
  left_join(prism_zip %>%
              select(-lat, -lng) %>%
              mutate(year = year - 1)) %>%
  group_by(region, desc, year) %>%
  summarise(
    yr_sol = mean(value),
    yr_temp = mean(tmean),
    yr_rain = mean(ppt),
    sd_sol = sd(value),
    sd_temp = sd(tmean),
    sd_rain = sd(ppt),
  ) %>% 
  filter(region %in% c(1, 2, 5)) %>% 
  ungroup() %>% 
  mutate(
    desc = str_remove(desc, "\\([[:alpha:]]+\\)")
  ) %>% 
  ggplot(aes(
    y = yr_temp,
    x = yr_sol,
    color = forcats::as_factor(region))) +
  geom_point() +
  geom_smooth(method = "lm") +
    scale_color_manual(
    values = c(
      "1" = "tomato2",
      "2" = "darkslategray4",
      "3" = "springgreen3",
      "4" = "brown",
      "5" = "goldenrod1",
      "6" = "gray50",
      "7" = "deeppink3",
      "8" = "gray17",
      "9" = "slateblue2"
    ),
    labels = c(
      "1" = "1: Desert",
      "2" = "2: Southeast",
      "3" = "3: High Plains",
      "4" = "4: Rainforest",
      "5" = "5: Arid Prairie",
      "6" = "6: Cold Desert",
      "7" = "7: Forested Mountains",
      "8" = "8: Fescue Belt",
      "9" = "9: Upper Midwest & Northeast"
    )
  ) +
  theme_classic() +
  facet_wrap(~ desc, 
             nrow = 3,
             scales = "free_y") +
    theme(
    plot.title = element_text(
      size = 22,
      face = "italic",
      margin = margin(
        t = 0,
        r = 0,
        b = 13,
        l = 0
      )
    ),
    axis.title = element_text(
      size = 16),
    axis.title.y = element_text(margin = margin(
      t = 0,
      r = 13,
      b = 0,
      l = 0
    )),
    axis.title.x = element_text(margin = margin(
      t = 13,
      r = 0,
      b = 0,
      l = 0
    )),
    axis.text = element_text(
      size = 16),
    legend.text = element_text(
      size = 14),
    strip.text.x = element_text(size = 14)
  ) +
  labs(x = "Mean CG solution for the year",
       color = NULL,
       y = "Mean temperature for the year",
       title = str_wrap("Mean annual temperature is not directly related to mean annual weaning weight CG solution",
                        width = 38))